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Abstract 

We consider the following type of reaction-diffusion systems, motivated by a specific 
problem in the area of heterogeneous catalysis: 

A pulse of reactant gas of species A is injected into a domain Pci' 1 that we refer to 
as the reactor. This domain is permeable to gas diffusion and chemically inert, except for 
a few relatively small regions, referred to as active sites. At the active sites, an irreversible 
first-order reaction A-* B can occur. Part of the boundary of U is designated as the reactor 
exit , out of which a mixture of reactant and product gases can be collected and analyzed 
for their composition. The rest of the boundary is chemically inert and impermeable to 
diffusion; we assume that instantaneous normal reflection occurs there. 

The central problem is then: Given the geometric parameters defining the configuration 
of the system, such as the shape of U, position of gas injection, location and shape of the 
active sites, and location of the reactor exit, find the (molar) fraction of product gas in the 
mixture after the reactor U is fully evacuated. Under certain assumptions, this fraction 
can be identified with the reaction probability —that is, with the probability that a single 
diffusing molecule of A reacts before leaving through the exit. 

More specifically, we are interested in how this reaction probability depends on the rate 
constant k of the reaction A B. After giving a stochastic formulation of the problem, 
we consider domains having the form of a network of thin tubes in which the active sites 
are located at some of the junctures. The main result of the paper is that in the limit, as 
the thin tubes approach curves in Mr, reaction probability converges to functions of the 
point of gas injection that can be described fairly explicitly in their dependence on the 
(appropriately rescaled) parameter k. Thus, we can use the simpler processes on metric 
graphs as model systems for more complicated behavior. We illustrate this with a number 
of analytic examples and one numerical example. 

Reaction-diffusion, metric graphs. 


1 Introduction 


Consider the system described schematically in Figure [0 The bounded domain U c rep¬ 
resents a chemical reactor packed with inert granular particles permeable to gas diffusion. At 
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the point x e U, one injects a small pulse of gas molecules of species A. The pulse diffuses in U 
with a certain diffusion coefficient D before escaping through the reactor exit, marked in the 
figure with a dashed line and denoted dU a . 

Now suppose a certain small subset A of U contains a catalytically active material promot¬ 
ing the irreversible first-order reaction A -* B with rate constant k. We call A the active zone 
Then a certain fraction of the molecules in the initial pulse may interact with A and convert 
to B before exiting the reactor. This number is called the fractional conversion. Motivated 
by certain problems in heterogeneous catalysis (see U® and references therein) we seek to 
determine the fractional conversion in terms of the rate constant k, the diffusion coefficient T> 
and the geometric configuration of A and U. 


reflecting boundary 



Figure 1: Schematic description of the chemical reactor. U is the domain of diffusion; the shaded 
regions represent active sites, whose union is denoted A. The boundary of U is the union of a reflecting 
part dU r and an absorbing part dU a - The point x is where reactant gas of species A is injected. We 
are interested in the probability that a gas molecule starting at x and undergoing Brownian motion will 
react before reaching dU a - 


The practical concern is that k may be difficult to determine experimentally, because the 
reaction A ->■ B may involve, at the microscopic level, a complex sequence of steps such as 
adsorption and surface reaction. On the other hand, measuring the fractional conversion is 
comparatively simple. Similarly, one can carefully control in a lab the size and shape of 
the reactor, the size and shape of the active zone, and the nature of the particulate matter 
supporting the diffusion (hence the diffusion coefficient T>). Determining k in terms of these 
other variables is therefore of considerable interest. 

We intend to set up a probabilistic model for determining k. Let us briefly describe what 
physical assumptions we impose for the model to work. First, we assume that the pressure 
within U is low enough, and the injected pulse is narrow enough, that the diffusion follows 
the Knudsen regime. In other words, the diffusion is driven at the microscopic level by tiny 
collisions between diffusing molecules and the particulate matter making up the reactor bed— 
not by collisions between molecules themselves. Second, we assume the pulse is small enough 
that any reactions occurring within U during the course of the experiment negligibly alter the 
composition of the active sites. These assumptions are physically reasonable and can in fact 
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be arranged for in a lab; see m or [10] for more information. 

With the set-up just described, write a(x ) for the fractional conversion, or ctk(x) when 
necessary to emphasize dependence on the rate constant k. Because of our assumption about 
the Knudsenian character of the diffusion, we can identify ak(x) with the reaction probability; 
that is, with the probability that a single molecule, entering the chamber at x and undergoing 
reflecting Brownian motion in U with diffusion coefficient T>, converts to B before hitting dU a - 
This identification is possible because every molecule’s path in the reactor is independent of 
every other’s. Therefore a single pulse of N molecules is equivalent to N independent pulses 
of one molecule. For this reason, we can model the reaction A -> B by killing a reflecting 
Brownian motion (RBM) in U with a rate function of the form r = kq where k is the reaction 
rate and q is the indicator function on A, or perhaps a smooth approximation thereof. 

By RBM in U, we mean a diffusion process (Xt, P x ) with sample paths in U and satisfying 
the stochastic differential equation 

X t - Xq = aB t + [\(X s )d£ s (1) 

Jo 

where a = V2T>, Bt is an ordinary d-dimensional Brownian motion, u is the inward-pointing 
unit normal vector field on dU, and £ s is a continuous increasing process (called the boundary 
local time) that increases only on {t : Xt e dU}. The SDE (JTJ) is called the Skhorohod or 
semimartingale decomposition of (X t ). See, for example, [5, [ 7 ] for more information about 
RBM’s that have a semimartingale decomposition. As for killing a diffusion, see [20]. 

The existence of an RBM (Xt) in U with a decomposition (JT]) requires that dU satisfy 
certain modest smoothness conditions. For example, C 2 smoothness is more than sufficient 
to ensure that (pQ) has a pathwise unique solution; see [H, HUH!, ESI El]• See [T6| for general 
notions about stochastic differential equations. Henceforth we shall assume dU is at least C 2 
unless otherwise stated. Furthermore, we assume the first hitting time of (Xt) to dU a is finite 
almost surely. As for A , we assume for now that it’s a compact subset of U with a regular 
boundary, and impose further conditions when necessary. 


It will be convenient to introduce the survival function i/>k(x), defined as the probability 
that Brownian motion starting at x reaches dU a without being killed (that is, without reacting). 
We also write if(x) when k is understood. Then the reaction probability we seek is ak(x) = 
1 - 'ipk( x )- Our model for the reaction and diffusion within U entails that 


ih(x) = 


exp 



kq(X s )ds 


(2) 


where T is the first hitting time of (Xt) to dU a , that is T := inf{t > 0 : Xt e dU a }- 


Functionals such as the right-hand side of © are well-known and have the following physical 
interpretation in terms of our catalysis problem: the probability that a gas molecule injected at 
x and following the sample path lo hasn’t reacted decreases exponentially with the amount of 
time spent in the chemically active region A, the reaction constant k being the exponential rate. 
Thus the expression exp j-A; / 0 T ^ q(X s (u)) dsj represents the probability, for that sample 
path, that the molecule does not react by the time it exits the reactor, and then ipk(x) is the 
expected value of this probability. 
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It is also well-known that ifk(x) can be expressed as the solution to an elliptic boundary 
value problem in U. Our numerical illustration below is based on this fact: 

Proposition 1. Suppose that u is a function which is bounded and continuous in U, satisfies 

DAu - kqu = 0 

in U \ dA, together with the boundary conditions ^ = 0 on dU r and u - 1 on dU a ■ Suppose 
furthermore that dA is a set of zero potential for RBM in U. Then u = if within U. 

Remark. Let S = inf{t > 0 : Xt € A} be the first hitting time of X to A. Then kq{X)ds 
is only positive on (S < T) = the event that X t hits A at all before leaving through dU a . 
Accordingly, letting k -*■ oo reduces © to ipoo(x) = P x [S < T] = the probability that, starting 
from x e U, the RBM X t escapes through dU a without ever hitting A. We write Ph(x) = 
1 - 'fioo(x) for the complementary probability that Xt hits A at all before leaving. This last 
probability can also be determined by solving a boundary value problem. Namely, if v(x) 
solves Av = 0 in V := U \ A, together with the boundary conditions u = 1 on dA, u = 0 on dU a , 
and n • S7v = 0 on dV \ [dU a u dA] then v(x) = Ph(x). 


Elsewhere we will address the existence of solutions satisfying the conditions in Proposition 
1. For now, we simply assume that a solution exists and can be approximated using the Finite 
Element (FE) method. This is not a serious defect since we only intend to use Proposition 1 
to corroborate numerically the predicted form of ccfc(x) in Theorem 2 below. 

To understand what to expect, then, consider the case where U is a line or a metric graph, 
and A = {a}, a single point. Then simple manipulations involving the strong Markov property 
(described in detail below) show that otk(x) factors as ak(x) = P^{x)f{x,k) where P\fix) is 
the hitting probability for A (described in the Remark under Proposition 1) and f(x,k ) is the 
reaction probability conditional on starting from a. Clearly P\,(x) doesn’t depend on k, so 
that f(k,x ) is the quantity of interest. In this case, f(x,k ) is simply equal to «fc(a), and it 
follows from basic properties of local times (described in detail below) that 


f(x,k ) 


A k 

1 + Xk 


(3) 


where A is a constant that doesn’t depend on k. For this reason we regard the formula <[3J) as 
“typical” behavior for situations where A can reasonably be thought of as a single point in the 
state space. This involves some coarse-graining, because single points are polar for Brownian 
motion. 


To illustrate this behavior, consider the domain in Figure [21 Taking the length of the 
shortest side of the polygonal region as unit, the radius of the disc-shaped active site is 0.1. 
Here, A is not a single point, but is nevertheless small enough relative to U that we can 
reasonably hope that a factorization of the form ctk(x) = Pt l (x)f(x,k) holds. That is, if we 
define f(x, k) by 


f(x,k) 


ak(x ) 

Ph(x) 


(4) 
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Figure 2: Example in dimension 2 of a reactor domain with a small active site. 


and then compute the right-hand members of © using the boundary value problems deter¬ 
mining atk(x) and Ph(x ), then we expect to see the behavior described in ©. 

To this end, we computed f(x, k ) as defined in © for a fixed x and 50 different equally- 
spaced values of k e [0,100] using FEniCS, and then “solved for A” under the working assump¬ 
tion that © holds. That is, we defined 

A , = ! f(x,k ) = 1 j> k (x)/P h (x) 

kl-f(x,k) k 1 - ip k (x)/P h (x) 

and checked A’s dependence on k. As expected, the resulting A’s exhibited little dependence 
on k . Denoting by A max and A min the maximum and minimum values of A over the range of 
fc’s we looked at, we found (A max - A min )/A min < 0.005 with mean value A = 0.02342. 

It will be shown that the above expression © is indeed a typical approximation of f(x, k) 
when the active region is a small single site. Roughly speaking, A is the amount of time 
accumulated at a, starting from a, before departing through the reactor exit. In particular, 
A depends only on the geometry of the reactor and the diffusion coefficient, and not k. For 
more general active site configurations, more complicated but related approximation formulas 
apply. 

We call the regions to which our main result will apply fat graphs. These are graph-like 
domains comprising a number of thin tubes and junctures that converge, in an appropriate 
sense, to an underlying metric graph T = (V, £). Here, £ denotes the edges of T, and V the 
vertices. Then each tube in the domain corresponds to an edge e e £, and each juncture 
corresponds to a vertex v e V. We assume that the active region A is concentrated in the 
juncture regions. The precise definitions are given in section [2} For now, figure [3] conveys a 
rough idea. 

There are two scale parameters involved in the approximation: e > 0, which gives the 
thickness of the tubes in the domain, and h > 0, which gives the radius of the active zone 
around each active site. As e j, 0, the domain collapses to the skeleton graph, and as h { 0, 
active sites collapse to vertices. At the same time, reaction activity, given by the rate constant 
k h = k/h must increase accordingly. 
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Figure 3: A fat graph in R d with its skeleton graph T. Edges have relative radii r e and e is a scaling 
parameter. Active sites, defined as regions where reaction can occur, have a relative radius S with 
scaling parameter h. The reaction constant must also be scaled as k h = k/h as h [ 0. 


Thus let f e,h (x,k/h ) denote the quantity introduced above that gives the dependence of 
the reaction probability on the rate constant. That is, 


f e,h (x, k/h) 


^k(x) 

Ph(x)' 


The main result is then the following. 


Theorem 2. Let U e be a bounded fat graph whose skeleton metric graph T has finitely many 
vertices and edges. Let C be the set of vertices corresponding to active sites of U e and 6 the 
relative radius of the active sites. Define n.= k5/T>, where D is the diffusion constant. Let |C| 
denote the number of active vertices. Then, under the hypotheses of Proposition [7] below, 

lim lim f e ’ h (x, k/h) = 1 - £ Pv( x )^j^, 
h^Oe^O ffi, A(/c) 


where p v (x) is the probability that a diffusing particle starting from x hits C for the first time 
at v, conditional on hitting C at all, and X(n), X v (k) are polynomials in n of degree at most 
|C| and |C| — 1, respectively. The coefficients of A(k),A„(k) depend only on geometric properties 
ofT: lengths of edges, degrees of vertices, location of exit and active vertices. When C = {u}, 
this limit reduces to 


lim lim f e ’ h (x, k/h ) 
h^>0 e ->0 


\n 

1 + Xk, 


It would be interesting to find methods that could tell a priori how the coefficients of A (k) 
and X v (k) can be expressed in terms of the geometry and topology of the graph, that is, as 
function of lengths, degrees, etc. A starting point would be to explore in a systematic way the 
properties of these polynomials for families of graphs. Here we are content with showing a few 
examples only. 


The reaction probability for the examples of figure |4] are easily obtained by solving the 
linear system of equations indicated in section [3 In all cases we assume that the point of gas 
injection is the left-most vertex. The results are as follows. 
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Figure 4: A few examples of graphs for which the reaction probability is easily computed. The vertex 
types are: • = inert, © = active, 0 = exit. The labels U indicate edge lengths and the Xi on the middle 
graph of the left column are the coordinates of the active nodes. On the third graph of the left column, 
n Q is the number of vertices of type 0. 


1. Top graph on the left column: 


o(k) 


2 Ik 


1 + ‘IIk 

Naturally, only the length of the edge between the active and exit vertices matters. On 
the other hand, if the initial edge is completely eliminated so that the starting point for 
Brownian motion is the active vertex, then the term 2 Ik should be replaced with Ik. 


2. Top graph on the right column: 

deg (v)Ik 

a(K) =-—— 

1 + deg (v)Ik 

where deg(u) is the degree of the active vertex. As in the first example, only the length 
of the edge connecting the active vertex to the exit matters, but the number of shorter 
edges leading to inert vertices influences a, regardless of those edges’ lengths, so long as 
their are greater than zero. 

3. Middle graph of the left column. Let lj = Xj - Xj-\. The value of o(k) can be obtained 

recursively as follows. Set g\ = 1 and gj+\ = gj + lj+i ( g\ h -f gj) k for j = 2, ... ,m. (The 

vertex 0 is at position x m+ i.) Then 

/ \ 9m+ 1 ~ 1 

a(K) = -. 

9m +1 

It is interesting to note the following property concerning optimal arrangement of active 
nodes: When k is small, 

cx(k) = [(L - cci) + —h (L - x m )] k + higher order terms in k 

where L is the total length of the graph. So if k is small, one maximizes reaction 
probability by clustering all the active nodes near the entrance point. On the other 
hand, for large values of k, 


cx(k) = hlylm+id 71 + lower order terns in k. 
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One easily obtains that the coefficient of K m attains a maximum when the active vertices 
are equally spaced. 


4. Middle graph of the right column: 


a(ft) = a(oo) 


Xk 

1 + \k 


where a(oo) = Zi/(Z q + h) and A = 2 Z(Zq + h)/(lo + l\ + l). 


5. Bottom graph of left column: 

2(Z o + n 0 Z)«: 

a(K) =-—-. 

1 + 2(Zq + HqI^k 

Despite having multiple active vertices, this system behaves, in its dependence on n, like 
one with a single active vertex. 


6. Bottom graph of right column: 

. . 3Z(1 + Zi/c)k 

a(K) =- 7~/ -:—r—• 

1 + 3Z (1 + l\ k) k, 

As expected for a system with two active vertices, this is the quotient of second degree 
polynomials in k. 


2 Metric graphs and fat-graph domains 


We introduce here notation and terminology concerning metric graphs, together with an as¬ 
sociated class of domains in M. d we call fat graphs. Basically, a metric graph is an abstract 
graph realized as a collection of curvilinear segments; a fat graph is a neighborhood of T whose 
boundary satisfies some smoothness conditions. 

More precisely, let S = (V, £) be an abstract graph with vertex set V and edge set £. We 
assume the edges in £ are oriented, and denote by s,t : £ -> V the source and target vertex 
functions. For each edge e e £ let e be the inverse of e, which is the same edge given the 
opposite orientation. Thus s(e ) = t(e) and t(e) = s(e ). We assume that £ is closed under the 
inverse operation. We also assume that |£|,|V| < oo where j • | indicates the cardinality. 

Associate to each ueVa point in M d , still denoted v (so that we now think of V as a subset 
of the Euclidean space) and to each e e £ a smooth curve 7 e : [0, l e ] M fi , parametrized by 
arclength, such that 7 e ( 0 ) = s(e ) and 7 e (Z e ) = t(e). Thus, Z e is the length of the curvilinear 
segment E e := 7 e ([0,Z e ]). We assume that max e Z e < oo and write l e = |e| = \E e \. We also set 
7e(s) : =7e(|e| -s). 

The union T := Uee£ E e will be called a metric graph , denoted T. With slight notational 
abuse we also write T° for U e e£ where E° = 7 e ((0,Z e )). Thus, T° = T \ V. The word metric 
refers to the natural distance defined by minimizing a path between two points. This induced 
metric makes T into a separable metric space. 
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Each edge has a natural coordinate y e = 7 J 1 : E e -> I e . By means of this coordinate we 
can identify functions on E e with functions on [0, Z e ]. Similarly, we can identify a function 
/ : r -> M with a collection of functions on the various coordinate intervals [ 0 , Z e ] by defining 
/ e (s) = f('y e (s)) for s e [0, Z e ]. Thus we have an obvious way of checking that / is continuous: 
each / e must be continuous on (0 ,l e ) in the ordinary sense, and the extensions to [0, Z e ] must 
agree at the vertices, in the sense that f e = / e 2 (0) whenever t(e i) = s(e 2 ). The set of 
continuous functions on T is denoted C'(T). 


Similarly, we can define the derivative of / e C'(T) at a point x = 7 e (s) e E° by (/ 0 7 e) , (s), 
when this derivative exists in the ordinary sense. At a vertex v, we define the one-sided 
derivatives 


( D e f)(v) := lim 
s|0 


/(7c(-s)J - /(-’I 
s 


when the limit exists. Thus D e f(v) is the directional derivative of / at v, pointing into e. 
Clearly this definition only makes sense for e e £ such that s(e) = v. Note that, in general, we 
do not require that the directional derivatives at v all agree. 


We now define the fat graph U e as a union of certain e-tubes U* (one for each e e £) 
together with e-junctures J* (one for each v e V). Some care is required in formulating the 
definitions; however, Figure [3] should convey the right idea. The following conditions ensure 
that the construction works properly: 


1. If e,e' are any two edges such that s(e) = s(e') = v, then T e (u) * T e '(u). 

2. There exists a ro > 0 such that T'(s) = 0 for s € [0,ro] and all e e £. 

3. The curves comprised by T do not intersect each other and have no self-intersection. 

4. Each 7 e is C 3 . 


Now fix an e > 0. We begin by defining the e-tubes corresponding to the edges E e in T. 
For each e e £, let there be given a relative radius r e > 0; and then, for each e > 0 sufficiently 
small, let 17* be a tubular neighborhood of E° with cross-sectional radius r e e. According to 
pL2] . U* has the nearest point property with respect to E°, meaning each x e U e has a unique 
point TT e (x) € E° nearest to x, and the induced mapping ir e : U* -» E° is a C 2 submersion. 
Furthermore, the distance function 

x i-> d(x,E° e ) := inf{y e : j|x-yj|} 

is C 3 near E° (because y e is C 3 — see [12])) so that the unit normal vector field, given by 

rT(x) = Vd(x,E e )/\\Wd(x,E e )\\, x e dU*, 
is C 2 away from the ends of U*. 

Now, the full U *’’s may intersect near the vertices. For this reason, we use assumptions [2] 
and [3] from above to introduce constants Co > 0 and eo > 0 which depend only on T and have 
the following property: if each U* is shrunken by a length cqc at its ends, then U* n U* = 0 
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whenever ei,e 2 e £ are distinct and e $ eo- In summary, we have a family {Uf : e e £, e < eo} 
such that: (i) for fixed e ^ eo the 17 £, s are pairwise disjoint; (ii) (J e e£ OJ I 18 - 8 the nearest point 
property with respect to r°, and the natural projection 7r £ onto T° is O 2 ; (iii) U e e£ 0 £ I T° as 
ejO. 

Next, we define the juncture regions J £ . For each v e V let there be given a “template” 
region such that 

Jv U U U? 

e:s(e)=n 

is simply-connected and has a boundary smooth enough that the unit normal vector field is 
still C 2 away from the ends. In other words, J v lines up smoothly with its incident eo-tubes. It 
is clear that such a J v exists. It is also clear that we can shrink or enlarge J v so that J v n Uf° 
is a cylindrical wafer of height (ro -co)eo > 0 and radius r e e o, for each e e £. Then J £ is defined 
in an obvious way by scaling down homothetically: 

J £ := {iel‘ i :u + e~ 1 (x - v) e J v } . 



Figure 5: The e-juncture at v in part of a fat graph domain U e with skeleton T. The r e are the 
relative radii of the tube cross sections. Note that the part of the edge curve y e contained in is 
straight. 


Finally, with the e-tubes and e-junctures as above, we define the fat graph with fatness 
parameter e and relative radii r e as 


U e 


U Jv 

veV 


U Ul 

_ee£ 


From the construction, the unit normal vector field on dU e is C 2 , and U e i T as e j 0. It is also 
clear that we can extend the projections 7r e to be defined inside the J £ ’s in such a way that 
7r e : U e -> r is continuous, and 


sup sup ||vr £ (x) - v\\ = 0(e). 

vzV xtJfj 

Furthermore, we have that 

sup ||7r e (ec) - x|| = O(e) (5) 

X£U e 

uniformly in x e U e . 
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3 Limit of diffusions from the fat graph to its graph skeleton 


In this section we review some background material related to diffusions on fat-graph domains 
and their limits as the fat graphs shrink down to their underlying metric graphs. The main 
result quoted here comes from [T], which extends the earlier work m, with modifications 
required for our needs. 

Let T be a metric graph in M. d , and U e a family of fat graphs with skeleton T, fatness 
parameter e > 0, and relative radii r e (e e £) as defined in Section [2] Let a = (cW(x)) be a 
dxrf matrix-valued function on and b = (&*(x)) a vector field in M d . We assume that a is 
uniformly positive definite, and that a and b are both bounded and Lipschitz continuous in 
W l . Define the differential operator 

L = l'ZaV(x)D i D j + '£b i (x)D i , 

where Di is partial differentiation in X{. 

Writing a(x ) = a(x)a l (x), with a assumed positive-definite and Lipschitz continuous, con¬ 
sider the stochastic differential equation 

-XI = £<7 (XI) dW s + fjbixi) ds + J* IT (XI) dt s , (6) 

in which W is a d-dimensional Brownian motion, n is the inward pointing unit normal vector 
field on dU e and (i\) is a continuous increasing process adapted to the filtration of X and 
increasing only on the set [t '■ Xf € dU 6 }. Since n e is C 2 , the geometric conditions in [21] are 
clearly met. Therefore (J6]) admits a strong solution, in the sense of Pi- 

Let us describe the behavior of (Xf) as e | 0. First, define an operator on C(T) as 
follows. For each e e £, let Cl(E° e ) denote the space of functions on with two bounded and 
continuous derivatives; and then let L e act on / e C 2 (F1°) as 

(L e f)(x) = i ||a t (x)T e (x)|| 2 /"(a) + [(b(x),T e (x)) + (a(x)T' e (x),T e (x))] f e (s ), 

where x = J e (s), T'(x) = (T e o j e )' (s) and (•,•) is the ordinary Euclidean inner product. To 
“paste together” the L e ’s, we must specify what happens at the vertices. Thus, we define 
1>(Lr) to contain those functions / e C(T) for which: 

1. f e e Cl(E° e ) for each e e £; 

2. for each v € V, and e e £ such that v = s(e), the one-sided limits lim s ;o(^e/)(7e(s)) exist 
and have a common value, denoted (Lf)(v); 

3. for each v e V, 

E Pv(e)(D e f)(v) = 0 (7) 

e:s(e)=v 
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where the numbers p v (e) are defined by 


Pv(e) 


r 


d -1 
e 


X/e':s(e ')=v 


/yd 1 

e' 


Then, for / e D(Lr) as above, we define L-pf by 


(L r )(x) 


{(L e f e )(y e (x)) if xtE° 

jlim y ^ x (L e /)(y) if x e V. 


( 8 ) 


According to Theorem 3.1 in [14| . there is a diffusion process (Xt) on T generated by 
(Lr,D(Tp)). On the other hand, Theorem 4.2 in [I] says that (Xf) converges in distribution 
to ( Xt ) as e | 0 if Xq converges in distribution to a T-valued random variable. Since pathwise 
uniqueness holds for ([6]) under our assumptions, “in distribution” can be replaced by “almost 
sure” in the preceding sentence. In particular, this is true if the starting point Xq is a point 
x € T for all e > 0. 

Theorem 3. Let U e be a fat graph with skeleton T. Assume that Xq = x for all e > 0 and 
some ieT, Then X e converges in distribution to X with initial value Xq = x. If pathwise 
uniqueness holds in (j6|) , then X e converges a.s. to X. 

Remark. The coefficients p v (e) appearing above have the following probabilistic significance: 
if t$ is the first exit time of X t from a ball of radius 5 at v, then p v (e) = lirn^o PufXfra) e E e \. 


4 Diffusions with killing 


Let U e be a fat-graph domain with skeleton T. The terminology and assumptions of Sections 
[2] and [3] remain in force. 

We introduce a function r(x) which represents the rate of killing of a diffusion in U e or 
in r. This function is assumed to depend on a parameter h > 0 whose role will become clear 
later; roughly, r h will “collapse” to a vertex condition when h goes to 0. But for the moment 
we suppress h in order to simplify the notation. 

Thus, let r : M cZ -» [0, oo) be a non-negative measurable function, representing the chemical 
reaction rate. We are interested in processes obtained from X e on U e (resp. X on T) by 
killing using the rate function r. Loosely, killing with rate r means forming new processes 
Y e (resp. Y) which behave like X e (resp. X) until Jq r{X e s ) ds (resp. Jq r(X s ) ds) exceed an 
independent exponential random time; after this time, they are sent to a cemetery state A. In 
this situation, the Feynman-Kac formula says that Y e and Y have extended generators 


L r f = Lf - r(x)f and L r r f = L r f - r(x)f, 

respectively, where L and Lr are as defined in section [3] and the domains are the same. Fur¬ 
thermore, the semigroup of Y e acts on functions as 


Pt.f(x) = E, 


exp 


f\(Xl)ds}f(Xt) 
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with a similar statement holding for Y. In particular, Pfl(x) = E x [exp {-/ 0 * r (X £ s ) ds}] is the 
probability that Y £ is still in U £ , i.e. still “alive” at time t. With this in mind, we define the 
survival function of the process at x e U £ as 


f> e (x) := E x 


exp 



(9) 


Here, T £ is a random time which can be thought of as the time that X e is absorbed at the 
reactor, as we now explain. 



-> 

efO 


a 


Figure 6: A portion of the absorbing part of the boundary of U e and its limit in T. The process X e in 
U e is killed when it reaches the dotted line separating the regions indicated in light and dark shading. 
The limiting process X on T is killed when it reaches a. 


Let there be given a certain subset V a = {aj} c V of degree 1 vertices in T to be regarded 
as the reactor exit. Write (with abuse of notation) dU £ for the closure of the union of the 
corresponding juncture regions in U £ , and dUf for the rest of the boundary of U £ . We call dUf 
the absorbing part of the boundary and dUf the reflecting part. Then dU £ f V a as e | 0. (See 
Figure [6l) Let Tf be the first hitting time of X £ to dU £ , and T a the first hitting time of X 
to V a . For now, we assume that T £ and T a are finite a.s., and that T £ -+ T a.s. Under these 
circumstances, we have the following: 


Proposition 4. Suppose that Tf and T are a.s. finite and that T £ T a a.s. If x ef then 


if e ' h (x)^E x exp | - r h (X s )ds 


as e | 0, where T a is the hitting time to V a of the limiting process X. 


With some additional smoothness on r h , this conclusion can also be recast as a boundary 
value problem, which will be useful for computations. The following conclusion is standard. 
See, for example, [2]. 

Corollary 5. Suppose that u is a bounded and continuous function on U e that solves the 
equation Lu - r h u = 0 in U e , together with the boundary conditions 

— = 0 on dUf and u = 1 on dU £ . 
on 

Then u = il> e ’ h , with if £ ’ h defined as in ([9]), and 

if £ ’ h (x) - fj h (x) 

for every x e T, where is the solution to the corresponding ordinary differential equation on 
the graph T, namely: L^u - r h u = 0 with boundary condition u = 1 on V a . 
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5 Collapsing active zones towards vertices of T 


In the previous sections we described a conservative diffusion X = ( Xt ) on a metric graph T 
generated by an operator Lp acting on a domain characterized by the vertex condition dTJ). 
This X was obtained by collapsing the conservative diffusion X e on the fat graph U e down to 
T as e i 0. We also described a nonconservative process Y h = ( Y t h ) obtained by killing X using 
the rate function r h . The killed process Y h has generator Lp - r h acting on the same domain. 
The function r h represents the rate of chemical activity on the active zones. 

Now, we wish to collapse the active zones, i.e. the regions in T where r h is positive, to a 
collection of vertices, as h | 0. We have in mind that all the chemical activity is concentrated 
on the active vertices as h \. 0, while the killing rate is increased towards oo in such a way that 
the overall effect is the same in the limit. 

For concreteness, we assume that the operator Lp acts as DJp- on each edge E e , where 
T> is the diffusion coefficient; also, that the active regions are balls B(v,h5 ) centered at some 
of the vertices with radius h5, where 5 > 0 is a fixed number with the units of distance, and 
h > 0 is a dimensionless parameter. Thus, each B(v,h5 ) is a star-shaped neighborhood of v 
consisting of v and a union of segments of edges of length h5 for each edge issuing from v. The 
killing rate function r h (x ) will then be assumed to have the form 

r h (x) = £ y 1 B(v,hS)^ ( 10 ) 

veV n 

where k v is either k or 0 depending on whether that vertex is to be considered active or not. 
It is known (as explained in Remark 2.5 of m that, as h l 0, 


Ex 


exp 



E, 


exp 


£««ACT a ) , 

veV J . 


where k v = kyd/T) and i v (T a ) is the semimartingale local time (as defined in [13]) at v evaluated 
at the hitting time T a . 

Write C for the set of active vertices in V , and let £e(t) denote the local time accumulated at 
C up to time t. In other words, £e(t) is the sum of the i v {t) for which k v is not zero. Then the 
limit in the above expression becomes [exp {-K^e(T a )}] where k = k5/T>. Also, the killed 
processes Y h converge to a process Y which is obtained from X as follows: run X until 
exceeds an independent rate 1 exponential; after this time, send X to the cemetery state A. 
In other words, kill X using the local time £&{ t ) rather than the integral /q r h (X s ) ds. Then 
the new process Y is has a generator which is defined in the same way as Tp for functions on 
T°. However, the domain of Y’s generator is characterized by a different vertex condition, as 
explained in the following: 

Proposition 6. Let r e,h (x ) be positive bounded functions on U e with the property that 

limr e ’ ft '(x) = r h \x ) 

ej.0 
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for all x € T, where r h (x) is defined at m. Let if e,h (x) and if(x) be the survival functions 
associated with Y e,h and Y, respectively; that is, 


^ h (x) = E a 


exp 


■,f Q Ta r h ’ e (X e a )ds 


fi(x) = K x [exp {-/c£ e (T a )}] 


If the conditions of Proposition [7] are met, then 

i / j ( x ) = limlim if e,h (x) 
h ->0 e ->0 

where k = k6/T>, as defined above. Furthermore, the generator of the limiting process has a 
domain characterized by the vertex condition 

Y p v (e)(D e u)(v) = k v u(v) (11) 

e:s(e)=v 

where k v = k if v e 6 and n v = 0 if v e V \ C. 


Sketch of proof. The limit statement has already been discussed. As for the vertex condition, 
we will show that a function in the domain of Y’s generator must satisfy m Thus, let 
( L,D(L )) be the generator of Y. Fix an active vertex » f G. Let Us = B(v,6 ) and fs the 
first exit time of Y± from Us- Also, write E„ for the law of Y± started at v. If F e D(L) then 
Dynkin’s formula (see [20] Ch. 3) reads: 


(LF)(v) = lim 

U 


E v [F(Y(t s ))]-F(v) 

E„[fs] 


Let e be an independent rate 1 exponential (we can always enlarge the probability space to 
accommodate such a variable) and then f = inf{t > 0 : n£e(t) > e}. Thus, f is the time that Y t 
jumps from T to the cemetery state A. Evidently f v = ts a £ where ts is the first exit time of 
X t from Us- In the denominator, then, we have 


~E„[t 5 ] = — 


as 5 i 0. On the other hand, Y t leaves Us either through a point ( 5 , e) (using local coordinates) 
if ts < C or by a jump to A if £ ^ t$. Since F( A) = 0 in the latter case, we have in the 
numerator, in view of the Remark under Theorem (3[ 


E v [F(Y(r s ))] - F(v) = Y, Fe(S)F v [T 6 < C] - F(v) 

e 


= Y F e(5)p v (e)K[T5 < C] -F(v) +o(6) 


Now (ts < C) = (n£ e (Ts) ^ e) = (k£ c (ts) ^ e) P„-a.s., because starting from v, only the local 
time at v can contribute to £e(t) before ts- Also, £ v (rs ) has an exponential distribution under 
P„, as explained below; we can compute its E„-mean as E„ [r^] = 5 using the methods in 
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Section [3 Therefore P v [k£ v (t§) ^ e] is the probability that a mean k5 exponential is less than 
an independent mean 1 exponential; this is easily found to be Therefore the numerator 

equals 

E Fe($)p v (e) - F(v) + o(5 ) - Y,Pv(e) [F e (5 - F(v) - k6F(v )] 

= 1 1 E^( e ) l D eF( v ) - kF (v)] 5 + o(5) 

1 + no g 

Existence and finiteness of the limit in Dynkin’s formula therefore requires that Y,e(D e F)(v) - 
kF(v) = 0 if F e D(L). Similar reasoning away from the vertex shows that L has to act as L 
there. For the rest of the details see HZ], □ 

Remark. The upshot of these considerations is that, in the limit as h [ 0 and then e | 0, we 
obtain an expression E x [exp{-K4(T a )}] for the survival probability. Therefore the problem 
of finding explicit formulas for reaction probabilities reduces to the problem of evaluating E x - 
moments of 4(T a ). There is a method for dealing with this issue in considerable generality, 
called Kac’s moment formula, which we explain in the next section. 


6 Explicit formulas for reaction probability 

We now focus on the dependence on n of the conversion probability a K (x ) = 1 - 'il } k( x ) for the 
reaction-diffusion process on a metric graph T, where the active sites consist of a set of vertices 
C c V. In this case, it is possible to obtain reasonably explicit formulas for a K (x) by simple 
arguments using the strong Markov property. First we consider the case where C is a single 
vertex: 

Proposition 7. Suppose C = {c} and T a = inf {t > 0 : X t € V a }. Set A := E c [£ c (T a )] = the 
expected local time accumulated at c up to time T , starting from c. Then 

\k 

a K (x) = Ooo(x)--—. (12) 

1 + AK, 

Proof. Write T c := inf{f : Xt = c}, the first time that Xt hits the active vertex. Since £ c (T a ) = 0 
on the event (T c ^T a ), we have 

a K (x) = E* [(1 -exp{-K4(T a )}) 1 (Tc< t q) ] 

= a 00 (x) -E x [exp{-«4(T a )} 1(t c <t 0 )] 

where a 00 {x) = E. x [l/ Tc<Ta )] is the probability that X t hits c before V a . Also, since £ c (t) does 
not start increasing until T c , we have that £ c (T a ) = £ c (T a ) o 0 Tc on (T c < T a ). Using the strong 
Markov property, and then the fact that X{T C ) = c a.s., we hnd that the right-most term in 
the second line of (fT3|) equals: 

Ex [exp {-k£ c (T a ) o 0 Tc } l (Tc<Ta) ] = E s [E x(Tc) [exp{-Av4(T a )}]l (Tc<Ta) ] 

= E c [exp{-«^ c (T a )}]E x [l (Tc<To )] 
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We arrive at 


a K (x) = a x {x) (1 -E c [exp{-«4(T a )}]). (14) 

Now, we use the fact that i c (T a ) must have an exponential distribution under P c . This follows 
from a simple argument which can be found on p. 106 of [T9]. Writing A := E c [4(4)] for the 
mean, and then explicitly computing the Laplace transform, we obtain the result 

E c [exp {~rJ: c {T a )}] = 1 

1 + Xk 

Inserting this last expression into (fl4l) establishes (fT2|) . □ 


When G consists of more than one point, the survival function takes the form 


A(x) = Ej, 



where 4(4) is the local time at the vertex c up to the exit time T a . In this case, it is still 
possible to obtain an explicit formula for if K (x) (hence a K (x )) because Kac’s moment formula 
can be used to evaluate the expectation appearing in fi K (x). 

Specifically, we can use the following corollary of Kac’s moment formula, which can be 
found in Section 6 of m 

Proposition 8. Write G = { Cj }*? 1 with C = ffG and define the Green’s function by 

Gij := g(ci,Cj) = E Ci [4 j(T a )]. 

Let Ki = n(ci) be a positive function on G. Then the function 


fi K (ci) := E Ci 


C 


ex P{- E KjZcjiTa)} 

3 = 1 


is the unique solution to the system of equations 


C 

1p(Ci) = 1 - E GijKjlf(Cj). 
3 =i 


In other words, if M K is the diagonal matrix with the entries of Kj along the main diagonal, 
then 

ii> K = (I + GM K )- 1 l e . (15) 

Remark. The expression (1151) for the survival function on Q can be recast as a formula involving 
determinants by using Cramer’s rule. To wit, let G^ denote the matrix obtained from G by 
subtracting row j from every row. (So, G^ has a row of 0’s in the j-th row.) Then 


; , , det (I + G&M k ) 
’ det(/ + GM«) 

See [19] Ch. 2 for additional details about (fT6l) . 


( 16 ) 
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Now we can repeat the same argument that was given in Proposition 0 replacing the single 
point c with the set C, so that T c becomes Tg = inf{i > 0 : Xt e C}. Then (11311 takes the form 

a K (x) = aoo(x) - E x [exp {-n£ e (T a )} o d Te l(T e<Ta )] 

where we abbreviate Xj=i ^c, (T 0 ) as te(T a ). To deal with the expectation in the last display, 
split the event (Tg < T a ) as U^i (X(Tq) = Cj,Tp < T a ), and then write 

Pj(x) = P x [X(T e ) = Cj\T e < T a ]. 

Thus, Pj(x ) is the probability of starting from x and first hitting C at Cj, conditional on hitting 
C at all. We get 


C 

E x [exp {-K£ e (T a )} o 0 Te t( Te <T a )\ = E F xi x (T e ) 


3 = 1 


= Cj, T e < T a ]E Cj [exp {-Ki e (T a )] 


C 


= Y,Pj( x )^x[Te < T a ]E Cj [exp{-^ e (r a )}] 
3 = 1 


Since P x [Tg < T a ] = aoo(x), we can apply (fl6l) with n(cj) all being the same constant k to 
express this last equation in the simple form 


a K (x) = aoo(x) 


!- T,Pj( x ) 


3 = 1 


det (I + kGW) 
det (/ + kG) 


(17) 


In the numerator, the determinant is a polynomial with degree ^ C - 1, because k appears in 
only C— 1 rows of the matrix. On the other hand the denominator is a polynomial with degree 
^ C. Furthermore the coefficients depend only on the G%j, which in turn depend only on the 
geometric properties of T. This explains the claim made in Theorem [21 


7 Remark concerning the examples 

In this section we present two slightly different approaches to computing the reaction probabil¬ 
ities. In both subsections, we let X = ( Xt ) be a diffusion process on T as described in section 
[21 For simplicity, we’ll assume that X is a Brownian motion with coefficient D, so that the 
operator Lp acts as Lpf(x) = T>-^f e (x) for x e E°. This X is a conservative diffusion process 
(i.e. without killing) which corresponds to the limit of a conservative diffusion process on U e 
as e 1 0. 

7.1 Method 1 

Previously, we explained how Kac’s moment formula can be used to express the reaction 
probabilities in terms of certain polynomials involving the coefficients Gjj = K Ci [£ Cj (T a )]. What 
remains, then, is to compute the Gij ’s from the graph and diffusion coefficients. To this end 
we apply the graph stochastic calculus from m- 
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Then we have, for each v e V, a process £ v (t) which is adapted to the filtration of X, 
continuous, increasing, and increases only on {t : Xt e v}. As explained in [13], t v (t) can be 
recovered from X by the formula 

where B(v,S ) is a ball around v of radius 6. Also, let C%(F) be the set of / e C'(T) with two 
continuous derivatives in T°. (The derivatives don’t have to extend to be continuous at the 
vertices.) Then m gives the following version of the Ito-Tanaka formula for T : 

Theorem 9. Let F e C^(T). Define, for each v e V, 

Pv( F ) = E Pv(e)(D e F)(v). 
ee£:s(e)=v 


Then 

F(X t ) - F(X 0 ) = M t+ f\L r F)(X s )ds + E Pv(F)£ v (t) 

J0 vtV 


where Mt is a continuous local martingale. 


Ito’s formula is all we need to find the G,j. Namely, to find E Ci [^ c .(T a )], we must find a 
function F e C%(T) for which (i) L^F = 0 on T°, (ii) F(v) = 0 for v e V a , and (iii) p v (F ) = 0 
except for v = Ci, Cj. In this case, (i) means that F is affine on each edge, so that F e (y) = a e y+h e 
for certain constants a e and h e and 0 ^ y ^ |e|, where we abbreviate \E e \ = |e| for the length 
of the segment E e . Choosing an orientation arbitrarily for each e, this gives a total of 2 E 
constants, with E = ffH, which must be determined compatibly with (i)-(iii). This reduces 
the problem to a discrete one on the underlying combinatorial graph. This method is best 
illustrated by working a few examples. In the examples, when it is not necessary to give names 
to the edges, we write [ui,U 2 ] for the edge with s(e ) = v\ and f(e) = v^. 

Example 1. Consider the graph shown in the upper right of Figure |4j Write c for the active 
vertex in the middle, a for the absorbing vertex at the far end, and v\,...,v n for the remaining 
inert vertices. For simplicity, we assume that the p c ([c, Uj])’s are all equal to -, i.e. that the 
relative radii were all equal. Starting from any Vj, the probability of hitting c before a is 1. So 
the conversion will have the form 


a, 


X v j) 


Xn 


1 + Ak 

where A = E c [4(F n )]. Orient the various edges so that c is the origin. Then we seek a function 
F which is linear on each segment and continuous on T. The condition p Vj (F) = 0 forces F 
to be constant on [c, Vj ] and we can choose this constant to be l, the length of [c, o]. The 
condition F(a ) = 0 means that F is linear with slope -1 on [c,a]. In other words, p c (F ) = 
Taking expectations in Ito’s formula, 


0-Z = 0h--A => A = nl. 

n 

From this we obtain the formula given in item 2 in the list of examples after Theorem [2] 


OLniVj) = 


nln 
1 + nln 
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Example 2. In a similar spirit, consider the second figure in the left column of Figure [4j Let v 
be the inactive vertex at the left, and ci,..., c m the active vertices, and a the absorbing vertex 
at the right. Starting from v, X hits C with probability 1, and with probability 1 this first hit 
occurs at ci. Thus pi(v) = 1 and Pj(v) = 0 for the other Cj and formula (1171) simplifies to 


a K (v) = 1 - 


det(I + kG^) 
det (/ + kG) 


To simplify, we again assume that P[ Cj , Ck ]( e ) = 1/2 for adjacent Cj,Ck and p v ([v,c\]) = 1. 
Now, to compute Gij, first write lj for the length of the segment to the right of Cj. Then set 
Lj = Eilj h- Thus, Lfc is the sum of the lengths of all the segments to the right of c*,. Regarding 
T as a single straight line, take F to be the function which is constantly Lj on [0,Cj] and then 
decreases linearly with slope -1 on [cj,a], so that F(a) = 0 and p Cj (F) = -1/2. Using this F in 
Ito’s formula shows that: if Cj ^ Cj then Gij = 2Lj , and if Cj > Cj then Gjj = 2L*. In particular, 
the matrix G^ must be strictly triangular, so that the determinant in the numerator of a K (v) 
is 1. Therefore we can write 

det(J + nG) - 1 
* clct(/ + kG) 


7.2 Method 2 

A different approach is to work instead with the process Y = (Yt) obtained from X by sending 
it to A at £, the first time that n£e(t) exceeds an independent rate 1 exponential. This new 
process Y is a non-conservative diffusion whose generator still acts as Lp, but whose domain 
now consists of functions F e C 2 (T) satisfying these vertex conditions: 

Y, Pv(e)(D e F)(c) = kF(c) ceC 

e--s(e)=c 

£ Pv (e)(D e F)(v) = 0 v€Vs[QuV a ] 

e:s(e)=v 


It follows that if F e G 2 (T) satisfies LpF(x) = 0 in T°, F(a) = 1 for a e V a , together with the 
vertex conditions (1181) . then F(Yt ) is a martingale. By optional stopping, 

F(v) = K V [F(Y (T a ))] = Pi;[7a < C] 

which means that F(v) = il> K (v), i.e. the survival function evaluated at v. 

Again we assume that p v (e ) = l/deg(u) whenever v = s(e); equivalently, that all r e ’s are 
equal. Since F must be affine on each edge, it is determined by its values on V, and its 
derivative DF can be regarded as a function on £: 

DF{e) = F(t(e))-F(s(e)) 

\e\ 

Therefore we can couch the problem of determining £ as a kind of discrete boundary problem 
on the combinatorial graph 9 rather than on the metric graph T. For this purpose, regard the 
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vertices in V a as the boundary of 9 j and the remaining vertices as the interior of 9 • Then F is 
determined by the equations 


E M li? ( r ( e )) = [ d eg(^K+ E l e l 1 1 ( 19 ) 

e:s(e)= v \ e:s(e)=v J 

for interior vertices and f(v) = 1 for exit vertices, v e V a . Reaction probability for the examples 
given in the introduction are easily obtained by solving the above system of linear equations. 

Remark. Equation (1191) can be regarded as a kind of discrete Feynman-Kac equation involving 
the discrete Laplacian on 9- See [6] for more information. 
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